/*******************************************************************************
	- Generate figures for robustness checks
	- Notes: Before running this do file, Regression_Robustness.do needs to be 
			executed first.
*******************************************************************************/


********************************************************************************
* initial settings
********************************************************************************

version 14.2
clear 
drop _all
set more off

*----------------------------------

// set your own local path where you put the JAERE_Submission folder
if c(username) == "zhenx" {
	global LOCALPATH "C:/Users/zhenx/Dropbox/Research/Temperature and Athletes/temp_acclimatization/Code" 
}
global PATH "$LOCALPATH/JAERE_Submission"
global DATA "$PATH/Data"

********************************************************************************
* Robustness Checks
* Figure A4, A5: nonlinear relation by event type
* Figure A6, A7: nonlinear relation by ore-competition temperature
********************************************************************************

*--------------------------------------------------------------------
* list of robustness checks
*--------------------------------------------------------------------
/*
	1: exclude covariates
	2: replace dew point with relative humidity
	3: add team-venue FE
	4: alternative clustered SE 
	5: no weather interpolation 	
	6: two-week pre-exposure
	7: use maximum temperature
*/

*** choose the robustness check
	local robust_case 1 // 1 2 3 4 5 6 7

	if `robust_case' == 1 {
		global REG "$PATH/Results/Regressions/Robustness/Exclude Covariate"
		global FIGURE "$PATH/Results/Figures/Robustness/Exclude Covariate"
	}
	if `robust_case' == 2 {
		global REG "$PATH/Results/Regressions/Robustness/Relative Humidity"
		global FIGURE "$PATH/Results/Figures/Robustness/Relative Humidity"
	}
	if `robust_case' == 3 {
		global REG "$PATH/Results/Regressions/Robustness/Add Team-Venue FE"
		global FIGURE "$PATH/Results/Figures/Robustness/Add Team-Venue FE"
	}
	if `robust_case' == 4 {
		global REG "$PATH/Results/Regressions/Robustness/Alternative Clustered SE"
		global FIGURE "$PATH/Results/Figures/Robustness/Alternative Clustered SE"
	}
	if `robust_case' == 5 {
		global REG "$PATH/Results/Regressions/Robustness/No Weather interpolation"
		global FIGURE "$PATH/Results/Figures/Robustness/No Weather interpolation"
	}
	if `robust_case' == 6 {
		global REG "$PATH/Results/Regressions/Robustness/Two Week Pre-Exposure"
		global FIGURE "$PATH/Results/Figures/Robustness/Two Week Pre-Exposure"
	}
	if `robust_case' == 7 {
		global REG "$PATH/Results/Regressions/Robustness/Maximum Temperature"
		global FIGURE "$PATH/Results/Figures/Robustness/Maximum Temperature"
	}

	cd "$FIGURE"

	*--------------------------------------------------------------------
	* nonlinear relation between temperature and performance
	*--------------------------------------------------------------------
	if `robust_case' != 6 {

		*****************************coefficient estimates******************************

		import delimited "$REG/Temperature_Effect_by_Event_Type.txt", clear

		input 
			60 0 0 0
		end

		gen j = "coef" 	 if v1 != ""
		replace j = "se" if v1 == ""
		carryforward v1, replace

		drop in 1/1
		drop if strmatch(v1,"*Constant*") | strmatch(v1,"*Obser*") | strmatch(v1,"*R-squa*") | ///
				strmatch(v1,"*Adjusted*") | strmatch(v1,"*Numbe*") | strmatch(v1,"*Robust*") | ///
				strmatch(v1,"*p<0.01*") | strmatch(v1,"*travel*")
		drop if v2 == ""
		nrow

		ren coef j
		ren VARIABLES var
		replace var = subinstr(var,".Meet_temp_bin","",.)
		destring var, replace

		if `robust_case' != 7 {
			label define var_label 40 "<40" 45 "40-45" 50 "45-50" 55 "50-55" 60 "55-60" 65 "60-65" 70 "65-70" 75 "70-75" 80 ">75"
		}
		if `robust_case' == 7 {
			label define var_label 40 "<40" 45 "40-45" 50 "45-50" 55 "50-55" 60 "55-60" 65 "60-65" 70 "65-70" 75 "70-75" 80 "75-80" 85 "80-85" 90 "85-90" 95 ">95"
		}
		label values var var_label

		global varlist Sprint Strength Endurance
		foreach v in $varlist {
			replace `v' = subinstr(`v',"*","",.)
			replace `v' = subinstr(`v',"(","",.)
			replace `v' = subinstr(`v',")","",.)
			destring `v', replace
		}

		reshape wide $varlist, i(var) j(j) string
		foreach v in $varlist {
			gen `v'H = `v'coef + 1.96 * `v'se
			gen `v'L = `v'coef - 1.96 * `v'se
		}

		grstyle init
		grstyle set plain, nogrid compact

		foreach v in $varlist {
			twoway ///
				(connected `v'coef var, lcolor(navy) lwidth(medthick) msize(medium) mcolor(navy) msymbol(circle)) ///
				(line `v'H var, lcolor(gs10) lwidth(medthick) lpattern(dash) cmissing(n)) ///
				(line `v'L var, lcolor(gs10) lwidth(medthick) lpattern(dash) cmissing(n)) ///
				, ///
				xlabel(40(5)80, valuelabel) ///
				ylabel(-0.2(0.05)0.05) ///
				yline(0, lpattern(solid) lcolor(red)) ///
				legend(size(small) order(1 "Estimate" 2 "95% CI") position(6) ring(0) region(fcolor(none) lcolor(none))) ///
				title(`v') ///
				xtitle("") ytitle("Performance Change Relative to 55-60F") ///
				xsize(6) ysize(6) 
			graph save "Temp_`v'.gph", replace
		}

		***************************temperature distribution***************************
		use "$DATA/TF_Data.dta", clear
		drop if Team_tavg <= 40
		drop if Meet_tavg == . | Team_tavg_7 == . | Meet_preciptotal ==. | Meet_ozone == .
		if `robust_case' != 7 {
			if `robust_case' != 5 {
				keep event_type Meet_tavg
			}
			if `robust_case' ==5 {
				keep event_type MeetNI_tavg
				ren MeetNI_tavg Meet_tavg
			}
			foreach event in Sprint Strength Endurance {
				histogram Meet_tavg if Meet_tavg >= 40 & Meet_tavg <= 75 & event_type == "`event'" ///
					, bin(45) fcolor(gs12) lcolor(gs14) xlabel(40(5)75) xsize(6) fysize(25) xtitle(Mean Daily Temperature (F)) ///
					saving("Meet_Tavg_Dist_`event'.gph", replace)
			}
		}
		if `robust_case' == 7 {
			keep event_type Meet_tmax
			foreach event in Sprint Strength Endurance {
				histogram Meet_tmax if Meet_tmax >= 40 & Meet_tmax <= 95 & event_type == "`event'" ///
					, bin(45) fcolor(gs12) lcolor(gs14) xlabel(40(5)95) xsize(6) fysize(25) xtitle(Max Daily Temperature (F)) ///
					saving("Meet_Tavg_Dist_`event'.gph", replace)
			}
		}

		********************************combined graph********************************

		foreach v in Sprint Strength Endurance {
			graph combine ///
				"Temp_`v'.gph" "Meet_Tavg_Dist_`v'.gph" ///
				, col(1) imargin(0 0.6 0 0) xsize(10) scale(1.2) 
			graph save "Temperature_Effect_by_Event_`v'.gph", replace
		}

		graph combine ///
			"Temperature_Effect_by_Event_Endurance.gph" ///
			"Temperature_Effect_by_Event_Sprint.gph" ///
			"Temperature_Effect_by_Event_Strength.gph" ///
			, col(3) imargin(0 0.6 0 0) xsize(10) scale(1.2)
		graph export "Temperature_Effect_by_Event.eps", replace

		foreach v in Sprint Strength Endurance {
			erase "Temp_`v'.gph"
			erase "Meet_Tavg_Dist_`v'.gph"
			erase "Temperature_Effect_by_Event_`v'.gph"
		}

	}

	*--------------------------------------------------------------------
	* nonlinear relation between temperature and performance 
	*--------------------------------------------------------------------

	*** generate a dataset to plot the figure using estimation results
		clear
		set obs 1000
		gen n = _n
		if `robust_case' != 7 {
			gen Meet_tavg = 30 + n*(90-30)/500 if n <= 500
			replace Meet_tavg = 30 + (n-500)*(90-30)/500 if n > 500
			keep if Meet_tavg >=40 & Meet_tavg <= 80
		}
		if `robust_case' == 7 {
			gen Meet_tavg = 30 + n*(95-30)/500 if n <= 500
			replace Meet_tavg = 30 + (n-500)*(95-30)/500 if n > 500
			keep if Meet_tavg >=40 & Meet_tavg <= 95
		}
		gen hot = cond(n<=500,0,1)

	*** set other control variables at the mean level
		gen wind = -0.0000445
		gen Meet_preciptotal = 0.0819
	    gen Meet_ozone = 0.0367131
	    gen travel_dist = 277.2341

	*** create splines
		if `robust_case' != 7 {
		    mkspline2 meet_cubspl= Meet_tavg, cubic di knots(49.6 56.7 63.8)
			mkspline meet_linspl1 50 meet_linspl2 60 meet_linspl3  = Meet_tavg, di
		}
		if `robust_case' == 7 {
			mkspline2 meet_cubspl= Meet_tavg, cubic di knots(59 67 75)
			mkspline meet_linspl1 60 meet_linspl2 70 meet_linspl3  = Meet_tavg, di
		}

	*** plot the figure 
		foreach model in Quadratic Cubic Cubic_Spline Linear_Spline {
			preserve
				estimates use "$REG/Endurance_`model'"
				predict yhat, xb
				predict sehat, stdp
				replace yhat = yhat 
				gen yhat_H = yhat + 2*sehat 
				gen yhat_L = yhat - 2*sehat 

				
				grstyle init
				grstyle set plain, nogrid
				twoway ///
					 (line yhat_H Meet_tavg if hot == 1, lwidth(medium) lcolor(red) lpattern(dash) cmissing(y) sort) ///
					 (line yhat_L Meet_tavg if hot == 1, lwidth(medium) lcolor(red) lpattern(dash) cmissing(y) sort) ///
					 (line yhat   Meet_tavg if hot == 1, lwidth(medium) lcolor(red) cmissing(y) sort) ///
					 (line yhat_H Meet_tavg if hot == 0, lwidth(medium) lcolor(blue) lpattern(dash) cmissing(y) sort) ///
					 (line yhat_L Meet_tavg if hot == 0, lwidth(medium) lcolor(blue) lpattern(dash) cmissing(y) sort) ///
					 (line yhat   Meet_tavg if hot == 0, lwidth(medium) lcolor(blue) cmissing(y) sort) ///
					 , ///
					 legend(order(1 "+/- 2SE" 3 "Warmest Pre-Exposure" 4 "+/- 2SE" 6 "Coolest Pre-Exposure") position(6) ring(0) region(fcolor(none) lcolor(none))) ///
					 xtitle("Temperature (F)") ytitle("Standardized Performance Relative to World Record") ///
					 xlabel(40(5)80) ///
					 ylabel(-2.6(0.1)-2.3) /// 
					 title(`model', size(medium)) ///
					 plotregion(margin(zero)) ///
					 scale(1.2) 
				graph save  "Endurance_`model'.gph", replace
			restore
		}

	*** combine the figures together
		grstyle set plain, nogrid 
		grc1leg2 "Endurance_Quadratic.gph" "Endurance_Cubic.gph" "Endurance_Linear_Spline.gph" "Endurance_Cubic_Spline.gph" ///
			, legendfrom(Endurance_Quadratic.gph) col(4) xsize(10) ysize(16)
		gr_edit .style.editstyle margin(zero) editcopy
		gr_edit .plotregion1.graph1.style.editstyle margin(medium) editcopy
		gr_edit .plotregion1.graph2.style.editstyle margin(medium) editcopy
		gr_edit .plotregion1.graph3.style.editstyle margin(medium) editcopy
		gr_edit .plotregion1.graph4.style.editstyle margin(medium) editcopy
		gr_edit .legend.style.editstyle margin(vmedsmall) editcopy
		gr_edit .style.editstyle declared_xsize(12.5) editcopy
		gr_edit .plotregion1.graph3.title.text = {}
		gr_edit .plotregion1.graph3.title.text.Arrpush Linear Spline
		gr_edit .plotregion1.graph4.title.text = {}
		gr_edit .plotregion1.graph4.title.text.Arrpush Cubic Spline
		graph export "Adapt_Function_Forms.eps", replace
		
	*** delete useless intermediate files
		foreach file in "Endurance_Quadratic.gph" "Endurance_Cubic.gph" "Endurance_Linear_Spline.gph" "Endurance_Cubic_Spline.gph" {
			erase `file'
		}




********************************************************************************
* Figure A8: nonlinear relation between temperature and performance of other events
********************************************************************************

global REG "$PATH/Results/Regressions/Sprint Strength"
global FIGURE "$PATH/Results/Figures/Sprint Strength"
cd "$FIGURE"

foreach event in Sprint Strength {

	*** generate a dataset to plot the figure using estimation results
		clear
		set obs 1000
		gen n = _n
		gen Meet_tavg = 30 + n*(90-30)/500 if n <= 500
		replace Meet_tavg = 30 + (n-500)*(90-30)/500 if n > 500
		gen hot = cond(n<=500,0,1)
		keep if Meet_tavg >=40 & Meet_tavg <= 80

	*** set other control variables at the mean level
		if "`event'" == "Sprint" {
			gen wind = 0.3642563
			gen Meet_preciptotal = .0803285
		    gen Meet_ozone = .0369447
		    gen travel_dist = 273.7413
		    global ylabel1 = -2.75
		    global ylabel2 = -2.5
		}
		if "`event'" == "Strength" {
			gen wind = .1348044
			gen Meet_preciptotal = .0794931
		    gen Meet_ozone = .0369213
		    gen travel_dist = 262.35
		    global ylabel1 = -4.6
		    global ylabel2 = -4.4
		}

	*** create splines
	    mkspline2 meet_cubspl= Meet_tavg, cubic di knots(49.6 56.7 63.8)
		mkspline meet_linspl1 50 meet_linspl2 60 meet_linspl3  = Meet_tavg, di

	*** plot the figure 
		foreach model in Quadratic Cubic Cubic_Spline Linear_Spline {
			preserve
				estimates use "$REG/`event'_`model'"
				predict yhat, xb
				predict sehat, stdp
				replace yhat = yhat 
				gen yhat_H = yhat + 2*sehat 
				gen yhat_L = yhat - 2*sehat 

				
				grstyle init
				grstyle set plain, nogrid
				twoway ///
					 (line yhat_H Meet_tavg if hot == 1, lwidth(medium) lcolor(red) lpattern(dash) cmissing(y) sort) ///
					 (line yhat_L Meet_tavg if hot == 1, lwidth(medium) lcolor(red) lpattern(dash) cmissing(y) sort) ///
					 (line yhat   Meet_tavg if hot == 1, lwidth(medium) lcolor(red)  cmissing(y) sort) ///
					 (line yhat_H Meet_tavg if hot == 0, lwidth(medium) lcolor(blue) lpattern(dash) cmissing(y) sort) ///
					 (line yhat_L Meet_tavg if hot == 0, lwidth(medium) lcolor(blue) lpattern(dash) cmissing(y) sort) ///
					 (line yhat   Meet_tavg if hot == 0, lwidth(medium) lcolor(blue)  cmissing(y) sort) ///
					 , ///
					 legend(order(1 "+/- 2SE" 3 "Warmest Pre-Exposure" 4 "+/- 2SE" 6 "Coolest Pre-Exposure") position(6) ring(0) region(fcolor(none) lcolor(none))) ///
					 xtitle("Temperature (F)") ytitle("Standardized Performance Relative to World Record") ///
					 xlabel(40(5)80) ///
					 ylabel($ylabel1(0.1)$ylabel2) ///
					 title(`model', size(medium)) ///
					 plotregion(margin(zero)) ///
					 scale(1.2) 
				graph save  "`event'_`model'.gph", replace
			restore
		}

	*** combine the figures together
		grstyle set plain, nogrid 
		grc1leg2 "`event'_Quadratic.gph" "`event'_Cubic.gph" "`event'_Linear_Spline.gph" "`event'_Cubic_Spline.gph" ///
			, legendfrom(`event'_Quadratic.gph) col(4) xsize(10) ysize(16)
		gr_edit .style.editstyle margin(zero) editcopy
		gr_edit .plotregion1.graph1.style.editstyle margin(medium) editcopy
		gr_edit .plotregion1.graph2.style.editstyle margin(medium) editcopy
		gr_edit .plotregion1.graph3.style.editstyle margin(medium) editcopy
		gr_edit .plotregion1.graph4.style.editstyle margin(medium) editcopy
		gr_edit .legend.style.editstyle margin(vmedsmall) editcopy
		gr_edit .style.editstyle declared_xsize(12.5) editcopy
		gr_edit .plotregion1.graph3.title.text = {}
		gr_edit .plotregion1.graph3.title.text.Arrpush Linear Spline
		gr_edit .plotregion1.graph4.title.text = {}
		gr_edit .plotregion1.graph4.title.text.Arrpush Cubic Spline
		graph export "Adapt_Function_Forms_`event'.eps", replace

}
	
*** delete useless intermediate files
	foreach model in Quadratic Cubic Linear_Spline Cubic_Spline {
		erase "Sprint_`model'.gph"
		erase "Strength_`model'.gph"
	}



